//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __statics_hpp__
#define __statics_hpp__

#include <vector>
#include <string>

#include "math3d.hpp"
#include "utils.hpp"
#include "screen.hpp"

class Screen;
class CoilData;

/** Generic static element that generates fields/occupies space. */
class StaticElement
{
public:
  virtual ~StaticElement ();

  /** Add fiel values at given position. */
  virtual void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield) = 0;

  /** Add the electrostatic potential from the given element. */
  virtual void addPotential (const vect3d& pos, prec_t& potential) = 0;

  /** Draw the elemnt onto a screen. */
  virtual void draw(Screen& screen, const Screen::color& c = Screen::BLACK) = 0;
  /** Get povray script of the given element. */
  virtual std::string get_povray_script () = 0;
  /** Get gnuplot script of the given element. */
  virtual std::string get_gnuplot_script () { return std::string(); }
  /** Test whether given position is in contact with the given current element. */
  virtual bool inContact (const vect3d& pos) = 0;
};

/** Integrated form of a current element where path is divided into tiny bits for calculations. */
class StaticIntegrated : public StaticElement
{
public:
  /** @internal */
  struct POSF
  {
    virtual void get_pos (prec_t t, vect3d& out) = 0;
    virtual void get_slope (prec_t t, vect3d& out) = 0;
    virtual ~POSF() { }
  };

  /** Set charge of the given integrate element for efield.
    * @param qpm coulombs per meter. */
  void setCharge (prec_t qpm);

protected:
  StaticIntegrated (prec_t currentAMPS, prec_t wire_radius, POSF *posf);
  virtual ~StaticIntegrated ();

  // StaticElement interface
  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  bool inContact (const vect3d& pos);

  void addPotential (const vect3d& pos, prec_t& potential) { }

protected:
  POSF *posf;
  prec_t currentAMPS;
  prec_t wr;

private:
  prec_t qpm;

  prec_t KmI;

  struct entry
  {
    vect3d middle;
    vect3d slope;
    prec_t qmod;
  };

  std::vector<entry> elist;

  void *calcmap_mem;
  entry *calcmap_begin;
  entry *calcmap_end;

  array16<prec_t> sse2_data;
  size_t sse2_consume;

  vect3d mid;
  prec_t gr;
  prec_t max_rm;

  vect3d plane_normal;
  prec_t pi_plane_d;
  prec_t min_d, max_d;
};

/** Current element formed by a cubic bezier curve. */
class StaticBezier : public StaticIntegrated
{
public:
  StaticBezier (prec_t currentAMPS, prec_t wire_radius,
      const vect3d& p0, const vect3d& p1, const vect3d& p2, const vect3d& p3);
  virtual ~StaticBezier ();

  std::string get_povray_script ();

private:
  vect3d p0, p1, p2, p3;
};

/** Base class for circular current elements. */
class RingBase : public StaticIntegrated::POSF
{
public:
  RingBase (prec_t currentAMPS, const vect3d& pos, prec_t radius, prec_t wire_radius, const vect3d& normal);
  virtual ~RingBase () { }

  std::string get_povray_script ();
  std::string get_gnuplot_script ();
  bool inContact (const vect3d& pos);
  void ringbase_draw (Screen& screen, const Screen::color& c);

protected:
  prec_t current;
  vect3d pos;
  prec_t radius;
  prec_t wire_radius;
  vect3d normal;
  prec_t plane_d;

  prec_t wrwr;
  prec_t min_rr, max_rr;
  vect3d nr;
  vect3d runit;
  vect3d lr;

  void get_pos (prec_t t, vect3d& out);
  void get_slope (prec_t t, vect3d& out);
  void get_closest (prec_t cosphi, vect3d& out1, vect3d& out2);

protected:
  vect3d p0;
  vect3d p0_pc;
  vect3d nXp0_pc;
};

/** Circle of current (or circle of charge for efield) with slow numeric integration. */
class StaticRingIntegrated : public RingBase, public StaticIntegrated
{
public:
  StaticRingIntegrated (prec_t currentAMPS, const vect3d& pos, prec_t radius, prec_t wire_radius, const vect3d& normal);
  virtual ~StaticRingIntegrated ();

  std::string get_povray_script ();
  bool inContact (const vect3d& pos);
};

/** Circle of current (or circle of charge for efield) with fast analytic calculation. */
class StaticRing : public RingBase, public StaticElement
{
public:
  StaticRing (prec_t currentAMPS, const vect3d& pos, prec_t radius, prec_t wire_radius, const vect3d& normal);
  virtual ~StaticRing ();

  // StaticElement interface
  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  std::string get_gnuplot_script ();
  bool inContact (const vect3d& pos);
  void setCharge (prec_t qpm);
  void addPotential (const vect3d& pos, prec_t& potential) { }

private:
  prec_t bfact;
  prec_t efact;
  prec_t charge;
  prec_t R, RR;
};

/** Base class for a static line. */
class StaticLineBase : public StaticElement
{
public:
  StaticLineBase (const vect3d& p0, const vect3d& p1, prec_t radius);

  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

protected:
  vect3d p0, p1;
  prec_t radius;
  vect3d l;
  prec_t l_len;
  vect3d ln;

  vect3d mid;
  vect3d lr;
  prec_t rr, dfms;
};

/** Line of current. */
class StaticLine : public StaticLineBase
{
public:
  StaticLine (prec_t currentAMPS, const vect3d& p0, const vect3d& p1, prec_t radius, prec_t rho = 0);

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);

private:
  prec_t currentAMPS;
  prec_t KmI;
  prec_t rho;
};

/** Infinite line of current. */
class StaticLineInfinite : public StaticElement
{
public:
  StaticLineInfinite (prec_t currentAMPS, const vect3d& p0, const vect3d& dir, prec_t radius);

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);
  void addPotential (const vect3d& pos, prec_t& potential) { }

private:
  prec_t currentAMPS;
  vect3d p0;
  vect3d normal;
  prec_t radius;
  vect3d lnormal;
  prec_t KmI2;
  prec_t rr;
};

/** Rectangular coil formed out of four line elements. */
class StaticRectangle : public StaticElement
{
public:
  StaticRectangle(prec_t currentAMPS, const vect3d& pos, const vect3d& right, const vect3d& up, prec_t wr);
  virtual ~StaticRectangle ();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);
  void addPotential (const vect3d& pos, prec_t& potential) { }

private:
  vect3d p1, p2, p3, p4;
  prec_t wr;
  StaticLine *l1, *l2, *l3, *l4;

  std::string make_povsegment (const vect3d& n1, prec_t d1, const vect3d& n2, prec_t d2,
      const vect3d& p1, const vect3d& p2);
};

class StaticCIBase : public StaticElement
{
public:
  StaticCIBase (const CoilData& cd, const vect3d& pos, const vect3d& normal, prec_t realI, prec_t realqpm);
  virtual ~StaticCIBase ();

  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);
  void addPotential (const vect3d& pos, prec_t& potential) { }

protected:
  const CoilData& data;
  vect3d pos, normal;
  prec_t realI, realQ;
  vect3d lr;
};

/** Coil interpolator, nearest neighbour. */
class StaticCINN : public StaticCIBase
{
public:
  StaticCINN (const CoilData& cd, const vect3d& pos, const vect3d& normal, prec_t realI, prec_t realqpm);
  virtual ~StaticCINN();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
};

/** Coil interpolator, nearest neighbour. */
class StaticCICU : public StaticCIBase
{
public:
  StaticCICU (const CoilData& cd, const vect3d& pos, const vect3d& normal, prec_t realI, prec_t realqpm);
  virtual ~StaticCICU();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
};

/** Coil interpolator, bilinear interpolation scheme. */
class StaticCIL : public StaticCIBase
{
public:
  StaticCIL (const CoilData& cd, const vect3d& pos, const vect3d& normal, prec_t realI, prec_t realqpm);
  virtual ~StaticCIL();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
};

class StaticLineCharge : public StaticLineBase
{
public:
  StaticLineCharge (const vect3d& p0, const vect3d& p1, prec_t rho0, prec_t rho1, prec_t blotradius = 0.0);

  void addPotential (const vect3d& pos, prec_t& potential);
  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);

private:
  prec_t rho0, rho1;
  prec_t rhodiff;
  prec_t rhodiff_by_llen;
  prec_t mod;
};

/** A charged sphere. */
class StaticChargedSphere : public StaticElement
{
public:
  StaticChargedSphere (const vect3d& pos, prec_t radius, prec_t q);
  virtual ~StaticChargedSphere ();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

private:
  vect3d pos;
  prec_t radius;
  prec_t q;
  prec_t qmod;
  prec_t radiusradius;
};

/** A sphere of uniform charge (charge within sphere). */
class StaticUniformSphere : public StaticElement
{
public:
  StaticUniformSphere (const vect3d& pos, prec_t radius, prec_t q);
  virtual ~StaticUniformSphere ();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

private:
  vect3d pos;
  prec_t radius;
  prec_t q;
  prec_t qmod;
  prec_t imod;
  prec_t R2;
  prec_t rho;
};

class Octree;

/** Octree based cube polywell interpolator. */
class StaticOctree : public StaticElement
{
public:
  StaticOctree (Octree& tree, bool do_bfield, bool do_efield, bool do_pot, prec_t Bmul = 1, prec_t Emul = 1);
  ~StaticOctree();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

protected:
  Octree &tree;
  bool do_bfield, do_efield, do_pot;
  prec_t Bmul, Emul;
};

class StaticOctree24 : public StaticOctree
{
public:
  StaticOctree24 (Octree& tree, bool do_bfield, bool do_efield, bool do_pot, prec_t Bmul = 1, prec_t Emul = 1);
  virtual ~StaticOctree24 ();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
};

class StaticOctree48 : public StaticOctree
{
public:
  StaticOctree48 (Octree& tree, bool do_bfield, bool do_efield, bool do_pot, prec_t Bmul = 1, prec_t Emul = 1);
  virtual ~StaticOctree48 ();

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
};

/** Static elements. */
class Statics
{
public:
  Statics ();
  ~Statics ();

  /** Add a new static element into the system.
    * Statics gets ownership of the given element object. */
  void addStaticElement (StaticElement *addElement);

  /** Get fields from given position. */
  void getFields (const vect3d& pos, vect3d& bfield, vect3d& efield);

  /** Test whether given position is in contact with any static elements. */
  bool inContact (const vect3d& pos);

  /** Set default ambiend bfield. */
  void setAmbientBField (const vect3d& abfield) { ambient_bfield = abfield; }
  /** Set default ambiend efield. */
  void setAmbientEField (const vect3d& aefield) { ambient_efield = aefield; }

  /** Draw all static elements onto given screen. */
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);

  /** Get povray script for all the static elements. */
  std::string get_povray_script ();
  /** Get gnuplot script for all the static elements. */
  std::string get_gnuplot_script ();

  /** Get integration step (meters). */
  static prec_t get_ddsq ();
  /** Set numerically calculate static elements integration step (meters). */
  static void set_ddsq (prec_t ddsq);

  /** Get the electrostatic potential. */
  prec_t getPotential(const vect3d& pos);

private:
  vect3d ambient_bfield, ambient_efield;
  typedef std::vector<StaticElement *> strips_t;
  strips_t strips;
};

#endif // !__statics_hpp__

